* Main estimation file for "Green Expectations: Current Effects of Anticipated Carbon Pricing".
* Runs primary regressions and saves results.


clear
capture log close
estimates clear

local dir="C:/yourdirectory"  /* main directory */
local covar_input = "`dir'"  /* directory with covariate files */
local contract_input = "`dir'"  /* directory with CL and dependent variable files */

log using "`dir'/logs/estimation.log", replace

set more off
set matsize 1000


*_______________________________________________________________________________


* Initialize

local covarlist reSP reLib reTnote reCR reBDI  /* covariates: don't include CL here because will vary by contract */
local covarfilelist SP Lib Tnote10 CR BDI  /* covariates: don't include CL here because will vary by contract */

local specif_list "base nocovar garch compositeindex"  /* specifications to run; "base" should always come first */
* possibilities: base, nocovar, garch compositeindex

local specif_list_tex "base nocovar garch"  /* specifications want to use in standard TeX documents */

local event_list "graham"  /* should read graham, but can set to "test" when want only one loop */

local autocorrel_check 1  /* =1: check for autocorrelation; =0: don't */


* These estimation window names are from old code; below we convert them into total trading days
local estwindow_list_graham "2 12" /* estimation windows, in months */ 

* Estimation windows in terms of trading days
local tradedays_2mo "60"
local tradedays_12mo "200"

local commod_list_graham_2mo "QL NG"  /* commodities */
local commod_list_graham_12mo "QL NG"  /* commodities */

local year_list_graham "2010 2011 2012" /* contract years to loop through for this event */

* Trading date of event day (e.g., Monday after weekend event) and surrounding days
local date_graham "04/26/2010"
local datebefore_graham "04/23/2010"
local dateafter_graham "04/27/2010"

* Define liquid contracts for each commodity in 2-month estimation window, making sure first two characters match commods variable
local QL2010liquid_2mo_graham "N Q U V X Z"
local QL2011liquid_2mo_graham "F G H J K M N Q U V X Z"
local QL2012liquid_2mo_graham ""
local NG2010liquid_2mo_graham "N Q U V X Z"
local NG2011liquid_2mo_graham "F G H J K M N Q U V X Z"
local NG2012liquid_2mo_graham "F G H J"

* Define liquid contracts for each commodity in 12-month estimation window, making sure first two characters match commods variable
local QL2010liquid_12mo_graham "V X Z"
local QL2011liquid_12mo_graham "F G H J K M N Q U V X Z"
local QL2012liquid_12mo_graham ""
local NG2010liquid_12mo_graham "V X Z"
local NG2011liquid_12mo_graham "F G H J K M N Q U V X Z"
local NG2012liquid_12mo_graham "F G H J"

* First contract will use
local first_year_graham "2010"
local first_month_2mo_graham "N"
local first_month_12mo_graham "V"

* Testing code
local estwindow_list_test "2" /* estimation windows, in months */ 
local commod_list_test_2mo "QL"  /* commodities */
local commod_list_test_12mo "QL"  /* commodities */
local year_list_test "2010" /* contract years to loop through for this event */
local date_test "04/26/2010"
local datebefore_test "04/23/2010"
local dateafter_test "04/27/2010"
local QL2010liquid_2mo_test "N Q U V X Z"
local QL2010liquid_12mo_test "V X Z"
local first_year_test "2010"
local first_month_2mo_test "N"
local first_month_12mo_test "V"

*_______________________________________________________________________________

* Begin looping

foreach event of local event_list { /* loop through events */

foreach estwindow of local estwindow_list_`event' { /* loop through estimation windows */

estimates clear

foreach commod of local commod_list_`event'_`estwindow'mo {  /* loop through commodities */

foreach specif of local specif_list{  /* loop through specifications */

use "`contract_input'/`commod'fut/`commod'`first_year_`event''`first_month_`estwindow'mo_`event''.dta", clear

local contract_tracker = 0  /* reset tracker */

foreach contract_year of local year_list_`event' { /* loop through contract years */

foreach contract_month of local `commod'`contract_year'liquid_`estwindow'mo_`event' {  /* loop through contract months */


local contract_tracker = `contract_tracker' + 1


* Get numeric year code for this contract month
local month_counter 1
foreach month_code in F G H J K M N Q U V X Z {
	if "`contract_month'"=="`month_code'" {
		local contract_number `month_counter'
	}
	local month_counter = `month_counter'+1 
}
if `contract_number' < 10 {  /* add 0 to front of single digits */
	local contract_number "0`contract_number'"
}
*display "`contract_number'"



if "`first_year_`event''`first_month_`estwindow'mo_`event''"~="`contract_year'`contract_month'" {  /* not first time through */
	quietly merge 1:1 time using "`contract_input'/`commod'fut/`commod'`contract_year'`contract_month'.dta", gen (d`commod'`contract_year'`contract_month')
	drop if d`commod'`contract_year'`j'~=3
	drop d`commod'`contract_year'`contract_month'
	
	* Add this contract's CL covariate
	if "`specif'" ~= "nocovar" {
		quietly merge 1:1 time using "`contract_input'/clfut/CL`contract_year'`contract_month'.dta", gen (dCL`contract_year'`contract_month')
		drop if dCL`contract_year'`contract_month'~=3
		drop dCL`contract_year'`contract_month'
	}
	
}
else {  /* first time through */
	
	* Define event window

	gen dummy_eventday = 1 if date == "`date_`event''"
	replace dummy_eventday = 0 if missing(dummy_eventday)
	
	gen dummy_daybefore = 1 if date == "`datebefore_`event''"
	replace dummy_daybefore = 0 if missing(dummy_daybefore)
		
	gen dummy_dayafter = 1 if date == "`dateafter_`event''"
	replace dummy_dayafter = 0 if missing(dummy_dayafter)

	* Add covariates
	if "`specif'" ~= "nocovar" {
		foreach covarfile of local covarfilelist {
			quietly merge 1:1 time using "`covar_input'/`covarfile'.dta", gen (d`covarfile')
			drop if d`covarfile'~=3
			drop d`covarfile'
		}
		quietly merge 1:1 time using "`contract_input'/clfut/CL`contract_year'`contract_month'.dta", gen (dCL`contract_year'`contract_month')
		drop if dCL`contract_year'`contract_month'~=3
		drop dCL`contract_year'`contract_month'
	}
	
	gen counter = _n
	gen estimation_window = 0 /* placeholder */
	gen days_from_event = 0 /* placeholder */
}


* generate counter for time variable; note that can differ for different covariates b/c e.g., LIBOR has British bank holidays
* returns have already been calculated in loaded files, though, using their own tsset data
sort time
replace counter = _n
tsset counter


/* Estimation window */	
gen target = counter if date == "`date_`event''"
egen event_day_count = min(target)
replace days_from_event = counter - event_day_count
replace estimation_window = 1 if abs(days_from_event)<=((`tradedays_`estwindow'mo'/2)+1) /* add 1 because event window includes day before and after */
replace estimation_window = 0 if estimation_window == .
*keep if estimation_window == 1  /* drop if outside of estimation window */
drop event_day_count target






if "`specif'" == "base" {
	
	gen unlog`commod'`contract_year'`contract_month' = exp(re`commod'`contract_year'`contract_month') - 1
	
	if `autocorrel_check' == 1 {
		quietly regress re`commod'`contract_year'`contract_month' dummy_eventday dummy_daybefore dummy_dayafter `covarlist' reCL`contract_year'`contract_month'  if estimation_window==1
		estat dwatson
	}
	
	ivreg2 re`commod'`contract_year'`contract_month' dummy_eventday dummy_daybefore dummy_dayafter `covarlist' reCL`contract_year'`contract_month'  if estimation_window==1, kernel(bartlett) bw(auto) small	
	
}
else if "`specif'" == "nocovar" {

	ivreg2 re`commod'`contract_year'`contract_month' dummy_eventday dummy_daybefore dummy_dayafter if estimation_window==1, kernel(bartlett) bw(auto) small

}
else if "`specif'" == "garch" {

	arch re`commod'`contract_year'`contract_month' dummy_eventday dummy_daybefore dummy_dayafter `covarlist' reCL`contract_year'`contract_month' if estimation_window==1, arch(1) garch(1) vce(robust) distribution(t) from(archb0 armab0) nonrtolerance ltolerance(1e-4) gradient gtolerance(999) /* trace */

}
else if "`specif'" == "compositeindex" {
	
	* Don't include dummies in this regression
	regress re`commod'`contract_year'`contract_month' `covarlist' reCL`contract_year'`contract_month'  if estimation_window==1

	* Save these for plotting and for composite index test
	/*
	predict resid_`commod'_`contract_year'`contract_month'  if estimation_window==1, residual
	drop if resid_`commod'_`contract_year'`contract_month' == .
	egen std_resid_`commod'_`contract_year'`contract_month'  = std(resid_`commod'_`contract_year'`contract_month' )
	*/
}


****** Store estimates *************


estimates store `specif'_`commod'_`contract_year'`contract_month'  /* will use in making tex file, do not need to specify event or est window in this label because will reset before then */


/* store residuals so can do SQ test and normality tests */
predict resid_`commod'_`contract_year'`contract_month' if estimation_window==1, residual
egen std_resid_`commod'_`contract_year'`contract_month'  = std(resid_`commod'_`contract_year'`contract_month' )


* store dummies' coefficients, standard erorrs, and p-values
if "`specif'" ~= "compositeindex" { /* have usual coefficients so store them */
	
	gen b_`commod'_`specif'_Day`contract_year'`contract_number' = _b[dummy_eventday]
	gen b_`commod'_`specif'_Before`contract_year'`contract_number' = _b[dummy_daybefore]
	gen b_`commod'_`specif'_After`contract_year'`contract_number' = _b[dummy_dayafter]

	gen se_`commod'_`specif'_Day`contract_year'`contract_number' = _se[dummy_eventday]
	gen se_`commod'_`specif'_Before`contract_year'`contract_number' = _se[dummy_daybefore]
	gen se_`commod'_`specif'_After`contract_year'`contract_number' = _se[dummy_dayafter]

	test dummy_eventday
	gen p_`commod'_`specif'_Day`contract_year'`contract_number' = r(p)
	test dummy_daybefore
	gen p_`commod'_`specif'_Before`contract_year'`contract_number' = r(p)
	test dummy_dayafter
	gen p_`commod'_`specif'_After`contract_year'`contract_number' = r(p)
	
}

} /* end of contract month looping */

} /* end of contract year looping */




/* save temp file b/c will want it for saving estimation results */
save "`dir'/temp_`specif'.dta", replace

/* save file with price and return time series for all of this commodity's contracts */
if "`specif'" == "base" {
	drop if estimation_window==0
	keep `commod'* unlog`commod'* date
	save "`dir'/output/`event'_timeseries_`commod'_`estwindow'mo.dta", replace
	use "`dir'/temp_`specif'.dta"  /* get all variables back */
}
/* save file with volume time series for all of this commodity's contracts */
if "`specif'" == "base" {
	drop if estimation_window==0
	keep vol`commod'* logvol`commod'* date
	save "`dir'/output/`event'_voltimeseries_`commod'_`estwindow'mo.dta", replace
	use "`dir'/temp_`specif'.dta"  /* get all variables back */
}



/* test for normality and store residuals for SQ test */
if "`specif'" ~= "compositeindex" {
	
	* drop missing values
	drop if estimation_window == 0
	
	* drop event window
	drop if dummy_eventday == 1
	drop if dummy_daybefore == 1
	drop if dummy_dayafter == 1	
	
	* Shapiro-Wilk test for normality
	swilk resid_`commod'_*	
	* saves in r(p)


	* Jarque-Bera test for normality
	foreach contract_year of local year_list_`event' { /* loop through contract years */
		foreach contract_month of local `commod'`contract_year'liquid_`estwindow'mo_`event' {  /* loop through contract months */		
			jb6 resid_`commod'_`contract_year'`contract_month'
			* doesn't save results?
		}
	}
	
	* save residuals for taking into Excel
	keep resid_*
	xpose, clear varname
	save "`dir'/output/`event'_resids_`commod'_`specif'_`estwindow'mo.dta", replace
	use "`dir'/temp_`specif'.dta"  /* get all variables back */	
	
}


/* Plot density of residuals */
if "`specif'" ~= "compositeindex" {
	
	* drop missing values
	drop if estimation_window == 0
	
	* drop event window
	drop if dummy_eventday == 1
	drop if dummy_daybefore == 1
	drop if dummy_dayafter == 1	
	
	* Kernel density
	foreach contract_year of local year_list_`event' { /* loop through contract years */
		foreach contract_month of local `commod'`contract_year'liquid_`estwindow'mo_`event' {  /* loop through contract months */		
			kdensity resid_`commod'_`contract_year'`contract_month', generate(x_resid_`commod'_`contract_year'`contract_month' y_resid_`commod'_`contract_year'`contract_month') nograph			
		}
	}
	
	keep x_* y_*
	save "`dir'/output/`event'_kdensity_`commod'_`specif'_`estwindow'mo.dta", replace
	
	use "`dir'/temp_`specif'.dta"  /* get all variables back */	
	
}


* Make matrix of coefficients, standard errors, and p-values for this commodity-specification pair
if "`specif'" ~= "compositeindex" { /* have usual coefficients so store them */

	keep if date == "`date_`event''"
	keep se_* b_* p_* date
	
	reshape long b_`commod'_`specif'_Day se_`commod'_`specif'_Day p_`commod'_`specif'_Day ///
		b_`commod'_`specif'_Before se_`commod'_`specif'_Before p_`commod'_`specif'_Before ///
		b_`commod'_`specif'_After se_`commod'_`specif'_After p_`commod'_`specif'_After, ///
		i(date) j(contract)

	save "`dir'/output/coeffs_`commod'_`specif'.dta", replace
}

}  /* end specification looping */


* Merge specification matrices for this commodity
foreach specif of local specif_list{ 
	if "`specif'" == "base" { /* are starting coefficient phase */
		use "`dir'/output/coeffs_`commod'_base.dta"
		erase "`dir'/output/coeffs_`commod'_base.dta"
	}
	else if "`specif'"~="compositeindex" {
		quietly merge m:m date using "`dir'/output/coeffs_`commod'_`specif'.dta"
		drop _merge
		save "`dir'/output/`event'_coeffs_`commod'_`estwindow'mo.dta", replace
		erase "`dir'/output/coeffs_`commod'_`specif'.dta"
	}
}


* Get merged dataset with estimation results
foreach specif of local specif_list{  /* loop through specifications */
	if "`specif'"=="base" {  /* first time through */
		use "`dir'/temp_`specif'.dta"
	}
	else {
		quietly merge 1:1 time using "`dir'/temp_`specif'.dta"
		drop _merge
	}		
}


* Make Tex table for each commodity's specifications
foreach contract_year of local year_list_`event' { /* loop through contract years */
foreach contract_month of local `commod'`contract_year'liquid_`estwindow'mo_`event' {  /* loop through contract months */
foreach specif of local specif_list_tex{  /* loop through specifications */

	estimates restore `specif'_`commod'_`contract_year'`contract_month'
	estimates store `specif'  											/* rename for column headings */
}
outreg2 [`specif_list_tex'] using "`dir'/output/tex/`event'/`estwindow'mo/`commod'_`contract_year'`contract_month'.tex", replace
}
}



* Composite index test for this commodity
foreach specif of local specif_list{  
	if "`specif'" == "compositeindex" {

		use "`dir'/temp_`specif'.dta", replace
		
		drop if estimation_window == 0

		tempname handle_compindex
		
		/***********************************************/
		* First do with unweighted average
		
		egen mean_resids = rowmean(std_resid_`commod'_*)
		
		* Regress on event dummies
		ivreg2 mean_resids dummy_eventday dummy_daybefore dummy_dayafter, kernel(bartlett) bw(auto) small

		
		file open `handle_compindex' using "`dir'/output/compindex_`event'_`commod'_`estwindow'mo.txt", text write replace
		
		file write `handle_compindex' "Simple mean " _n
		
		test dummy_eventday
		file write `handle_compindex' "Event day "
		file write `handle_compindex' "coefficient: " (_b[dummy_eventday]) ", stdrd error: " (_se[dummy_eventday]) ", p-value: " (r(p)) _n
		/* gen p_`commod'_`specif'_Day = r(p) */
		
		test dummy_daybefore
		file write `handle_compindex' "Day before event "
		file write `handle_compindex' "coefficient: " (_b[dummy_daybefore]) ", stdrd error: " (_se[dummy_daybefore]) ", p-value: " (r(p)) _n
		/* gen p_`commod'_`specif'_Before = r(p) */
		
		test dummy_dayafter
		file write `handle_compindex' "Day after event "
		file write `handle_compindex' "coefficient: " (_b[dummy_dayafter]) ", stdrd error: " (_se[dummy_dayafter]) ", p-value: " (r(p)) _n _n
		/* gen p_`commod'_`specif'_After = r(p) */
		
		
		
		
		
		/***********************************************/
		* Now do with weighted average (GLS)
		
		mkmat std_resid_`commod'_*, matrix(matrix_resids)
		matrix vector_ones = J(colsof(matrix_resids),1,1)
		matrix vector_composite_index = inv(vector_ones'*inv(matrix_resids'*matrix_resids)*vector_ones)*(vector_ones'*inv(matrix_resids'*matrix_resids)*matrix_resids')

		* Using svmat to convert the matrix into a variable
		matrix vector_composite_index = vector_composite_index'
		svmat vector_composite_index, names(composite_index) 
		
		* Regress composite index on event dummies
		ivreg2 composite_index dummy_eventday dummy_daybefore dummy_dayafter, kernel(bartlett) bw(auto) small

			
		file write `handle_compindex' "Weighted mean " _n
		
		test dummy_eventday
		file write `handle_compindex' "Event day "
		file write `handle_compindex' "coefficient: " (_b[dummy_eventday]) ", stdrd error: " (_se[dummy_eventday]) ", p-value: " (r(p)) _n
		/* gen p_`commod'_`specif'_Day = r(p) */
		
		test dummy_daybefore
		file write `handle_compindex' "Day before event "
		file write `handle_compindex' "coefficient: " (_b[dummy_daybefore]) ", stdrd error: " (_se[dummy_daybefore]) ", p-value: " (r(p)) _n
		/* gen p_`commod'_`specif'_Before = r(p) */
		
		test dummy_dayafter
		file write `handle_compindex' "Day after event "
		file write `handle_compindex' "coefficient: " (_b[dummy_dayafter]) ", stdrd error: " (_se[dummy_dayafter]) ", p-value: " (r(p)) _n
		/* gen p_`commod'_`specif'_After = r(p) */
				
		
		/*
		gen b_`commod'_`specif'_Day = _b[dummy_eventday]
		gen b_`commod'_`specif'_Before = _b[dummy_daybefore]
		gen b_`commod'_`specif'_After = _b[dummy_dayafter]

		gen se_`commod'_`specif'_Day = _se[dummy_eventday]
		gen se_`commod'_`specif'_Before = _se[dummy_daybefore]
		gen se_`commod'_`specif'_After = _se[dummy_dayafter]
		
		keep if date == "`date_`event''"
		keep se_* b_* p_* date time

		replace date = "`commod'"
		save "`dir'/output/compindex_`event'_`commod'_`estwindow'mo.dta", replace
		*/
				
		
		* Drop event window for normality test, SQ test, etc
		drop if dummy_eventday == 1
		drop if dummy_daybefore == 1
		drop if dummy_dayafter == 1	
				
		* Test for normality		
		swilk composite_index
		jb6 composite_index
		
		* save for kdensity (will overwrite in a second)
		save "`dir'/output/`event'_kdensity_`commod'_`specif'_`estwindow'mo.dta", replace
		
		* SQ test: save composite index and simple mean for Excel
		keep composite_index mean_resids
		xpose, clear varname
		save "`dir'/output/`event'_resids_`commod'_`specif'_`estwindow'mo.dta", replace
	
		* Kernel density		
		use "`dir'/output/`event'_kdensity_`commod'_`specif'_`estwindow'mo.dta", replace
		kdensity composite_index, generate(x_composite_index y_composite_index) nograph			
		keep x_* y_*
		save "`dir'/output/`event'_kdensity_`commod'_`specif'_`estwindow'mo.dta", replace
	
		file close `handle_compindex' 
		
	}
}



* Clean up
foreach specif of local specif_list{  /* loop through specifications */
	erase "`dir'/temp_`specif'.dta"
}


}  /* end of commodity looping */

} /* end of estimation window looping */

} /* end of event looping */



log close


